Bayesian Spatial Sampling Design

This paper presents a model-based approach to the problem of the optimal choice of a spatial sampling design in the presence of uncertainty about the distribution of the observations. Our approach, which is based on the polar spectral decomposition of stationary isotropic random fields, makes it possible to represent the spatial error process by a linear regression model with (infinitely many) uncorrelated random coefficients having expectation zero. Approximating this model by a finite number of terms we arrive at a random coefficient regression model for which well-known methods from convex analysis can be applied which allow the construction of efficient numerical (steepest descent) algorithms to find optimal spatial designs. We report on initial implementations of our spatial design algorithms aiming at the construction of spatial designs minimizing weighted versions of the mean squared prediction errors over a grid of prediction points. The ideas and implementations are illustrated for a real data set of Cs137 measurements from the region of Gomel, Belarus. INTRODUCTION Spatial sampling design, which means the “optimal” allocation of sampling points to spatial coordinates in order to improve spatial estimation and prediction in a well-defined sense, is a topic which has not found wide-spread attention so far in the spatial statistics literature. With the exception of the detailed account in Müller (2007) and a brief account in Diggle and Ribeiro (2007) this topic is not present in current geostatistics textbooks and monographs. The basic difficulty for rigorous model-based approaches to spatial sampling is the fact that the spatial correlation of the observations leads to analytically intractable objective functions. In this paper we will use a spectral representation for isotropic random fields which leads to an approximation by linear regression models with random coefficients. J. PILZ AND G. SPÖCK Contrary to the well-known Karhunen-Loeve-representation which has been used by Fedorov and Müller (1997) to tackle the spatial design problem, this spectral representaion is constructive and allows the effective computation of the approximating terms. Our approximation then allows the application of powerful tools from convex optimization theory and the well-developed experimental design theory for Bayes linear regression models to calculate spatial sampling designs for random fields. THE SPATIAL MIXED LINEAR MODEL We consider a mean square continuous (m.s.c.) random field {Z(x) : x ∈ D ⊆ R2} such that Z(x) = f (x)T β + ε(x), Eε(x) = 0 (1) where f (x) is a known vector of regression functions, β ∈ Rr a vector of unknown regression parameters and Cov(Z(x)),Z(y)) = C(x,y); x,y ∈ D (2) It is well-known that, under the assumption of compactness of D the error process of the m.s.c. random field has a uniformly convergent representation (in the m.s. sense) ε(x) = ∞ ∑ j=1 √ λ jα jψ j(x) for all x ∈ D, (3) where E(α j) = 0, Var(α j) = 1 and Cov(αi,α j) = 0 for i 6= j; i, j = 1,2, . . . The positive real numbers λ j stand for the eigenvalues and the continuous real functions are the corresponding eigenfunctions of C(·), which then has a uniform and (absolute) convergent representation C(x,y) = ∞ ∑ j=1 λ jψ j(x)ψ j(y) The components corresponding to the largest eigenvalues contribute most to the variability of the error process. Ordering the eigenvalues such that λ1 ≥ λ2 ≥ . . . and truncating the Karhunen-Loeve expansion at a sufficiently small eigenvalue λm we get an approximation of the random field (1) in the following form Z(x)≈ f (x)T β + m ∑ j=1 α jg j(x)+ ε0(x) (4) 2 GEOSTATS 2008, Santiago, Chile BAYESIAN SPATIAL SAMPLING DESIGN where g j(x) = √ λ jψ j(x),E(α j) = 0,Cov(αi,α j) = δi j; i, j = 1, . . . ,m and it is assumed that the residual process ∑j=m+1 √ λ jα jψ j(x) can be sufficiently well approximated by a white-noise process with some variance σ2 0 . The right-hand side of (4) is exactly in the form of a mixed linear model, with fixed parameters β = (β1, . . . ,βr) and random effects α = (α1, . . . ,αm) . This expansion of the covariance function was used by Fedorov (1996) and Fedorov and Flanagan (1997) in their work on experimental design. For practical applications, however, it will be very hard to solve the generalized eigenvalue problem which makes the KL-expansion not very attractive as a starting point for spatial sampling design. When restricting attention solely to trend estimation, i.e. estimation of EZ(x) = f (x)T β , then the results in Fedorov and Müller (2005) can be used for optimal design. It is, however, not obvious, how their approach can be applied to solve the more complicated design problem for prediction. POLAR SPECTRAL REPRESENTATION From now on we will additionally assume that the m.s.c. random field is covariance-stationary and isotropic, i.e. C(x,y) = C(||x− y||) for all x,y ∈ D. Then, according to Yaglom (1987), the covariance function can be represented in the form C(t) = ∫ ∞ 0 J0(tω)dS(ω), t ≥ 0 (5) where J0 is the Bessel function of the first kind and order 0, t stands for the length of the lag vector and S(·) is the so-called (polar) spectral distribution function associated with C(·). As such, S(·) is positive, monotonically increasing and bounded from above. On the other hand, knowing C(·) the spectral distribution function can be obtained from the inversion formula S(ω+)−S(ω−) 2 = ∫ ∞ 0 J1(tω)ωC(t)dt (6) where S(ω+) and S(ω−) denote the rightand left-hand side limits at ω and J1 denotes the Bessel function of first kind and order 1. Approximating S(·) by means of a step function with positive jumps δi = S(ωi+1)− S(ωi) at preselected points ωi, i = 0,1, . . . ,n− 1; where ω0 = 0, and changing to polar coordinates (t,φ) = (radius, angle), the polar spectral representation theorem for m.s.c. covariance-stationary and isotropic random fields tells us that the error process may be approximated as VIII International Geostatistics Congress 3 J. PILZ AND G. SPÖCK ε(t,φ)≈ ∞ ∑ m=0 { cos(mφ) n ∑ i=1 Jm(ωit)Umi } + ∞ ∑ m=1 { sin(mφ) n ∑ i=1 Jm(ωit)Vmi } where all the random variables Umi and Vmi are uncorrelated with variances Var(Umi) = Var(Vmi) = dmδi; and dm = 1 for m = 0, dm = 2 for m ≥ 1. Again, by truncating the above series at a sufficiently large term m = M, we get an approximation of our random field in form of a mixed linear model Z(x)≈ f (x)T β +g(x)T α + ε0(x) (7) The essential advantage over the previous KL-representation is, however, that the above representation allows the explicit computation of the additional regression function g(x) and the covariance matrix of the random vector α . Noting that Cov(ε(t1,φ1),ε(t2,φ2))≈ M ∑ m=0 dmcos(m(φ2−φ1)) n ∑ i=1 δiJm(t1ωi)Jm(t2ωi) it becomes clear that the components of the additional regression vector g(·) are made up of the following radial basis functions (cosine-sine-Bessel-harmonics) gm,i(t,φ) = cos(mφ)Jm(ωit); m = 0, . . . ,M; i = 1, . . . ,n gm,i(t,φ) = sin((m−M)φ)Jm−M(ωit); m = M +1, . . . ,2M; i = 1, . . . ,n BAYESIAN SPATIAL LINEAR MODEL Starting from our spatial mixed linear model (7) we may gain further flexibility with a Bayesian approach incorporating prior knowledge on the trend. To this we assume that the regression parameter vector β is random with E(β ) = μ ∈ Rr, Cov(β ) = φ (8) This is exactly in the spirit of Omre (1987) who introduced Bayesian kriging this way. He used physical process knowledge to arrive at “qualified guesses” for the first and second order moments, μ and φ . Some general ideas on specifying prior distributions for regression parameters may be found in Pilz (1991), ch. 3. On the other hand, the state of prior ignorance or non-informativity can be modelled by setting μ = 0 and letting φ−1 tend to the matrix of zeroes, thus passing the “Bayesian bridge” to universal kriging, see Omre and Halvorsen (1989). 4 GEOSTATS 2008, Santiago, Chile BAYESIAN SPATIAL SAMPLING DESIGN Now, combining (7) and (8), we arrive at the Bayesian spatial linear model (BSLM) Z(x) = h(x)T γ + ε0(x) (9) where h(x) = ( f (x) g(x) ) ,γ = ( β α ) ,Eγ = ( μ 0 ) =: γ0,Cov(γ) = ( φ 0 0 A ) =: Γ Here ε0(x) is white-noise with variance σ2 0 as before and A denotes the covariance matrix of α , resulting after the polar spectral approximation of the random field. Alternatively, g and A could also be chosen according to a truncated Karhunen-Loeve expansion. Since in the approximated model (9) we have uncorrelated errors, the prediction of Z(x0) at an unsampled location x0 ∈ D is equivalent to Bayes linear trend estimation, Ẑ(x0) = h(x0) γ̂B. Thus, collecting our observations at given locations x1, . . . ,xn ∈ D in the vector of observations Z = (Z(x1), . . . ,Z(xn)) and denoting by H = (hi(x j)) = (F ..G), where F = ( fi(x j)), G = (gi(x j)) are the usual design matrices formed with the regression functions f and g, respectively, we obtain Ẑ(x0) = h(x0) (HT H +σ2 0 Γ −1)−1(HT Z +σ2 0 Γ γ0) (10) Using the Sherman-Morrison-Woodbury matrix inversion formula the predictor (10) may be rewritten as Ẑ(x0) = f (x0) μ +(GAg(x0)+Fφ f (x0)) [σ2 0 In +GAG T +FφFT ](Z−Fμ) In this form, the analogy with the Bayesian kriging predictor ẐBK(x0) = f (x0) μ + (c0 + Fφ f (x0)) [K + FφFT ]−1(Z − Fμ) derived by Omre (1987) becomes obvious: GAg(x0) and σ2 0 In + GAG T approximate the vector of covariances c0 = (Cov(Z(x0),Z(x1)), . . . ,Cov(Z(x0),Z(xn))) and the covariance matrix K = Cov(Z) of observations, respectively. The total mean squared prediction error in the BSLM (9) is known to be E(Ẑ(x0)−Z(x0)) = σ(1+h(x0) [HT H +σ2 0 Γ]h(x0)). (11) Using the same matrix inversion formula as above, the right-hand side of (11) may be written as σ2 0 +g(x0) T Ag(x0)+ f (x0) φ f (x0) VIII International Geostatistics Congress 5 J. PILZ AND G. SPÖCK −(GAg(x0)+Fφ f (x0)) [σ2 0 In +GAGT +FφFT ](GAg(x0)+Fφ f (x0)), which approximates the Bayesian kriging error derived by Omre (1987): E(ẐBK(x0)−Z(x0)) = σ2 + f (x0) φ f (x0)− (c0 +Fφ f (x0)) [K +FφFT ](c0 +Fφ f (x0)) along the same lines of reasoning as above. Thus, we have derived the important result that trend estimation in the approximating BSLM (9) is equivalent to Bayesian kriging in the original model. SOME BASIC IDEAS OF REGRESSION DESIGN THEORY We are now going to discuss the problem of spatial sampling design. Very often, one is interested in minimizing some functional of the (Bayesian) kriging error, i.e. the total mean squared prediciton error (11) of Ẑ(x0) (which reduces to the universal kriging variance when letting φ−1→ 0). Here we consider the minimization of the weighted integral of this error over the area D⊂ R2 of interest, i.e. ∫ D E[Ẑ(x)−Z(x)]2w(x)dx→Min dn (12) where w(·) is a nonnegative function on D weighing the importance of the design points. The minimization is over all possible designs dn = (x1, . . . ,xn) ∈Dn of size n. Inspecting (11), this problem is equivalent to minimizing the functional ∫ D h(x)T [HT H +σ2 0 Γ −1]−1h(x)w(x)dx→Min dn (13) over D. This is, however, in the form of a standard Bayesian regression design problem, for which powerful tools and results are given in Pilz (1991). In (13), the design dn only enters in the design matrix H. The matrix M(dn) = 1 n H T H is the so-called information matrix of the design; for normally distributed observations it is proportional to the Fisher information matrix. Correspondingly, MB(dn) = 1 n [HT H +σ2 0 Γ −1] (14) is called the Bayesian information matrix. Denoting 6 GEOSTATS 2008, Santiago, Chile BAYESIAN SPATIAL SAMPLING DESIGN U = ∫ D h(x)h(x)T w(x)dx, (13) becomes equivalent to minimizing the trace functional Q(dn) := tr UMB(dn)→Min dn (15) The design dn = (x1, . . . ,xn) can be interpreted as a probability measure ξ giving weight n to each of the points xi ∈ D; i = 1, . . . ,n. Now, when embedding the exact designs dn = (x1, . . . ,xn) into the space of all discrete design measures ξ = {(x1, p1), . . . ,(xm, pm)} giving weight pi to the point xi ∈ D; i = 1, . . . ,m; m ≥ 1, and defining MB(ξ ) = m ∑ i=1 pih(xi)h(xi) + σ2 0 n Γ−1 (16) then the extended functional Q(ξ ) = tr[UMB(ξ )−1] has the nice property of being convex and continuous w.r.t. the convex set of all possible Bayesian information matrices gererated by some ξ . Therefore, powerful results from convex optimization theory can be used to compute optimal discrete design measures ξ ∗ which then can be approximated by some exact design d∗ n . ITERATION PROCEDURES We are now going to formulate iteration procedures for the construction of approximately optimal exact designs. Contrary to the construction of optimal discrete designs, which, among other design algorithms, are discussed in Pilz (1991), in this case we cannot prove convergence of the exact designs to the function value of an optimal exact design; we can only guarantee stepwise improvement of a given exact starting design, i.e. the sequence of function values decreases monotonically. We will give here two iterative algorithms, where the first serves for the construction of an initial design. Generation of an initial design The initial design is a one-point design which minimizes the design functional (15) among all designs of size n = 1. Note that such a design exists since the Bayesian information matrix is positive definite even for designs of size n = 1. Step 1. Choose x1 ∈ D such that Q(x1) = Min x∈D Q(x), where Q(x) = tr UMB(x) = tr U [h(x)h(x)T +σ2 0 Γ −1]−1 and set d1 = (x1). VIII International Geostatistics Congress 7 J. PILZ AND G. SPÖCK Step 2. Beginning with i = 1, find xi+1 such that Q(di +(xi+1)) = Min x∈D Q(di +(x)) and form di+1 = di +(xi+1). Continue with i replaced by i+1 until i+1 = n. Step 3. If i+1 = n then stop and take dn = (x1, . . . ,xn) as an initial design. Exchange Type Algorithm The second algorithm is an exchange type algorithm improving n-point designs and starting from an initial design. At each stage of the algorithm, the design point to be added is chosen such that it lies in the direction of the steepest descent of the design functional. Step 1. Use some initial design dn,1 = (x1,1, . . . ,xn,1) ∈ Dn of size n. Step 2. Beginning with s = 1, form the design dn+1,s = dn,s +(xn+1,s) by adding the point xn+1,s = arg in f x∈D Q(dn,s +(x)) to dn,s. Then form d j n,s = dn+1,s− x j,s by deleting that point x j,s from dn+1,s for which Q(dn+1,s− (x j,s)) = min 1≤i≤n+1 Q(dn+1,s− (xi,s)) Step 3. Repeat Step 2 until the point to be deleted is equivalent to the point to be added. ILLUSTRATION Now we are going to discuss an example demonstrating our approach for spatial sampling design. We use the so-called Gomel data set, a set of Cs137 measurements taken from a monitoring network in the region of Gomel (Belarus), an area in the neighbourhood of Chernobyl. We have logtransformed the Caesium137 concentrations at the given 591 data locations and calculated a weighted least squares fit to the empirical variogram of the logtransformed concentrations. Fig. 1 shows the corresponding covariance function and a worst case approximation to it based on the representation of the approximating covariance function. Here we have selected the largest angular frequency M = 35 and 68 radial frequencies ωi, i = 1, . . . ,68. The 2-dimensional spectral distribution function was approximated at these frequencies ωi by a step function and the stochastic process was approximated as well by Cosine-Sine-Bessel surface harmonics with random amplitudes. From the above frequency parameters M and ωi resulted a spatial Bayesian linear regression model with 4828 regression 8 GEOSTATS 2008, Santiago, Chile BAYESIAN SPATIAL SAMPLING DESIGN functions, which were used to approximate the error process. As variance for the error of the approximating Bayesian regression model we used σ2 0 = 0.3. For modeling the trend of the random function only a constant was used. Concerning the design problem, we were interested in minimizing the average total mean square error of prediction over a regular grid of 30× 30 points. This approximation was used, since the integrals for calculating the matrix U in the design problem tr(UMB(ξ )−1)→ min ξ , U = ∫ D ( 1 g(x) )( 1 g(x) )T w(x)dx here simplify to simple sums. We were interested in calculating 6-, 12-, 18-, 24-and 50-point designs, that should be added to the already available 591 sample locations. For the calculation of the desired exact designs we used a slight modification of the algorithms in the previous section. For that purpose, in every step of the algorithm two design points were added to the already available 591 sample locations by means of the ’initial design algorithm’. Thereupon, by means of the exchange type algorithm suboptimal 591+2-, 591+4-,...,591+50-point designs were generated step by step, such that in every step we got nearly optimal designs. Fig. 2 shows a suboptimal design of size 591+24. Obviously the algorithm selects design points in such areas at first, where the sampling density is low, a feature that we expect from good design algorithms. Moreover, every design point has been selected only once, a feature that seems to be a property of the polar spectral approximation, but is not a property of the used design algorithm. 0 50 100 150 200 250 300 −0.5 0 0.5 1 1.5 2 2.5 Figure 1: Covariance function and its approximation ACKNOWLEDGMENT This work was partially funded by the European Commission, under the Sixth Framework Programme, by the Contract N. 033811 with DG INFSO, Action Line IST-2005-2.5.12 ICT for Environmental Risk Management. The views expressed herein are those of the authors and are not necessarily those of the European Commission. VIII International Geostatistics Congress 9 J. PILZ AND G. SPÖCK −150 −100 −50 0 50 100 150 −150 −100 −50 0 50 100 150 24 Figure 2: Suboptimal 591+24 point design, Bayes Kriging 0 5 10 15 20 25 30 35 40 45 50 0.82 0.84 0.86 0.88 0.9 0.92 0.94 average kriging mean square error Figure 3: Decrease of the average squared prediction error when adding design points REFERENCESDiggle, P and Ribeiro, P (2007). Model-based Geostatistics. Springer, New York.Fedorov, V (1996). Design of spatial experiments: model fitting and prediction, vol. 13. Handbook ofStatistics, North-Holland.Fedorov, V and Flanagan, D (1997). Optimal monitoring network design based on mercer’s expansionof the covariance kernel. In Journal of Combinatorics, Information and System Sciences, ,no. 23.Fedorov, V and Mueller, W (2005). Optimum design for correlated processes via eigenfunctionexpansions. Tech. rep., to appear in J. Statist. Planning Inf.Mueller, W (2007). Collecting Spatial Data. Optimum Design of Experiments for Random Fields.Springer.Omre, H (1987). Bayesian kriging-merging observations and qualified guess in kriging. InMathematical Geology, vol. 19, pp. 25–39.Omre, H and Halvorsen, K (1989). The bayesian bridge between simple and universal kriging. In Math.Geology, vol. 21, no. 7, pp. 767–786.Pilz, J (1991). Bayesian Estimation and Experimental Design in Linear Regression Models. Wiley.Pukelsheim, F (1993). Optimal Design of Experiments. Wiley.Yaglom, A (1987). Correlation Theory of Stationary and Related Random Functions. Springer. 10GEOSTATS 2008, Santiago, Chile